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Abstract 

We present a numerical strategy to compute ensemble averages of 
coarse-grained two-dimensional membrane-like models. The approach 
consists in generalizing to these two-dimensional models a one-dimensional 
strategy exposed in [Blanc, Le Bris, LegoU, Patz, JNLS 2010, [6]], which 
is based on applying the ergodic theorem to Markov chains. This may be 
considered as a first step towards computing the constitutive law associ- 
ated to such models, in the thermodynamic limit. 

1 Introduction 

A standard problem in materials science is the computation of canonical av- 
erages. As an example, consider a system of P particles, at positions u = 
{u^, . . . , u^) 6 K^"^, and assume that its energy is given by a two-body interac- 
tion: 

Given an observable A : M:"^ — ^ R, and /3 > the inverse of the temperature, 
the canonical ensemble average of A is defined by (see e.g. [11]) 

A(w)e-'^^(")dM 
{A) = , (1.2) 
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where C is the macroscopic domain in which the particles are assumed to 
he. 

Such canonical averages relate the macroscopic properties of a system to 
the elementary phenomena at the microscopic scale, which are modelled by the 
atomistic energy For instance, the bulk pressure in a fluid is given by (|1.2p . 

for the observable 

1 ^ rIF' 

' ' 2— 1 

where T is the temperature and p the density. If the atoms are non-interacting, 
then E = Q and we recover Mariotte's law for ideal gas. In general, the rela- 
tion (|1.2p (with the above choice for A) allows to compute the pressure of a 
liquid at a given density and temperature, and hence to obtain its macroscopic 
constitutive law, based on the microscopic model (|l.ip . 

The integrals in (|1.2|) involve 3P variables, and this number is extremely 
large when one considers a macroscopic sample of matter. Therefore, a di- 
rect computation of {A) , by standard quadrature rules (using Gauss integration 
points), is not tractable. Several approaches have been proposed to deal with 
such problems (see e.g. f9" for a theoretical and numerical comparison of those). 
The approach we will follow here is to use the so-called overdamped Langevin 
dynamics. This stochastic dynamics reads 

du^ -VE{u)dt+ y^2/l3 dBt, (1.3) 

where Bt is a Brownian motion in dimension 3P. This dynamics is hence a 
steepest descent (according to the forces —'VE{u)), modified by a random noise, 
the magnitude of which is scaled by the temperature: the larger the temperature, 
the stronger the noise. The dynamics ()1.3|) is interesting because of the following 
ergodicity result: under general assumptions on the energy E, we have 



hm i r A[u{t))dt 

T->00 1 Jg 



(1.4) 



for (almost) all initial conditions u(0). In the long time limit, the random num- 

I r 

ber — / A{u{t))dt converges to a deterministic number, which is the canonical 

"'0 

average of interest. 

Thus, in order to compute (jl.2p , a standard approach consists in simulating 
the evolution of the system according to (11.31) (this can be done even for large 
systems), and averaging the quantity A{u{t)) along the obtained trajectory. The 
equation (|1.4p ensures that the time average of A{u{t)) converges to the desired 
quantity {A) in the long time limit. In practice, the equation (|1.3p is numerically 
integrated with the forward Euler scheme 



with the time step At, where Gm is a vector of Gaussian random variables in 
dimension 3P (all components of Gm are independent from each other, and they 
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are all distributed according to a Gaussian law of mean and variance 1). In 
turn, following (|1.4p . the canonical average (|1.2I) is approximated by 



1 *^ 



We described above how to compute the canonical average (jl.2[) using the 
overdamped Langevin dynamics (|1.3I) . We point out that we give no physical 
meaning to this dynamics, which is only considered here as a numerical tool to 
sample the Gibbs measure and allow efficient computation of canonical averages. 

As pointed out above, several such numerical tools have been proposed in the 
literature, and among them, the Langevin equation, which is closer to physical, 
Hamiltonian dynamics than the overdamped Langevin dynamics (jl.Sp . The 
Langevin dynamics reads 

du = M^^pdt, 

dp = -VE{u)dt-jM~^pdt+ y/¥iJl3 dBt, ^ ' 

where u G M.'^^ represents the position of the atoms, p G M'^^ their momentum, 
M is a mass matrix, and 7 > is a parameter. This dynamics shares the same 
ergodicity property as the overdamped Langevin dynamics: for any 7 > 0, the 
convergence (|1.4p holds along trajectories of ()1.5p . 

Note that, when 7 = 0, we recover from (|1.5|) the standard Hamiltonian 
dynamics (Newtonian dynamics). When 7 > 0, the Langevin dynamics includes 
some friction and some random noise, which is, as for the overdamped Langevin 
dynamics, scaled by the temperature. Finally, in the limit 7 00, and up to 
some time rescaling, the Langevin dynamics ()1.5|) converges to the overdamped 
Langevin dynamics (|1.3|) (see e.g. [24l Sec. 2.2.4]). 

In many cases, the observable A of interest does not depend on all positions 
M*, but only on a few of them. Our aim here is to derive simpler methods for 
computing canonical averages in such a case, based on a coarse-graining strategy. 
Such methods were used in [6] for one-dimensional models (see also |23J), and 
we adapt here the strategy to a two-dimensional membrane-like model. 

One possible approach to address such question is the QuasiContinuum 
Method (QCM). It was first introduced in [36l|37], and further developed in [211 
[28l [29l [33l [34l [38] , in the zero temperature case. It has been theoretically stud- 
ied in [11 El [3l 12 [3 [El [Ml [HI [III [Ml EB 135 See also [H [22] for some 
review articles. The method was adapted to the finite temperature case in [15] . 
resulting in a coarse-graining strategy to compute ensemble averages. Similar 
strategies were also proposed in [lOl [25] . 

In essence, any coarse-graining approach aims at using the fact that A de- 
pends only on the positions of some representative atoms at positions Ur G ri^'' , 
with Pr <C P. This is usually done by writing formally 

{A) = -^ f A(u,)e-'5^(")dM = ^ I A(M,)e-^^^'^("'-)du„ 
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where 



Zp 




exp {-f3E{u)) du, 




cxp {-^EcGiur)) du 



and where 



EcciUr) 




exp {—j3E{u)) diir 



(1.6) 



is the coarse-grained energy. In (|1.6p . the notation dur indicates that we inte- 
grate over all variables except Ur- The question then reduces to computing the 
coarse-grained energy Ecg, which can be considered as the free energy of the 
system, for the collective variable (or reaction coordinate) Ur- 

Note that Ecaiur) depends on P, and is in general difficult to compute for 
a given value of Pj. and P. However, one may hope that its expression simplifies 
in the limit when 



This is the case for instance in dimension one, in the case of nearest-neighbour 
interactions, as it is recalled in Section [Ol below. 

Rather than fixing P and letting P ^> oo, as described above, another asymp- 
totic regime is to keep P and Pr fixed and let fi go to -l-oo (vanishing temper- 
ature). This is the regime considered by the QCM approach, which is based, 
roughly speaking, on using a harmonic approximation of E{u) in (|1.6p . This 
then allows to compute explicitly Ecciur). 

The sequel of this article is organized as follows. In Section [O] we consider 
one-dimensional chains of atoms, and recall our results of [6]. We next present 
in Section [1.21 the two-dimensional membrane-like models we consider here, and 
the questions we are after. Section [2] is dedicated to the derivation of our 
numerical strategy, which is very much inspired by our previous works on one- 
dimensional models. The outcome of this derivation is an algorithm, that we 
have summarized in Section [2. 31 Numerical results obtained with this algorithm 
are presented (and compared with the reference model results) in Section [31 

1.1 Dimension one 

In [6^, we proposed a methodology to deal with one-dimensional models. Al- 
though the method allows to handle models with any finite range interaction, 
we recall it here in the special case of nearest-neighbour and next-to-nearest- 
neighbour interaction, refering to [6l [23] for more details. 

Consider a chain of P atoms at positions € R, with an energy of nearest- 
neighbour type: 



where h = 1/P is the average spacing of the atoms, and is assumed here to be 
equal to the characteristic length of the potential (as the equilibrium length of 



Pr is fixed, P — > cxD. 




(1.7) 
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z I—)- Wiz) is assumed to be of order 1, the equilibrium length of z i~> W(z/h) 
is of the order of h). To remove the translation invariance of the energy, we set 
uO = 0. 

As explained above, we consider the canonical average of observables de- 
pending only on a few atoms of the system. To simplify the setting, we assume 
that A depends only on . Thus, the quantity we wish to compute is 

p-i 



where Zp = / exp 



A (u^) exp 



,1+1 



du. 



p-i 



du. More precisely, we are going 



to compute the limit of {A)p when P ^ oo, that is in the thermodynamic limit, 
when the number of atoms present in the system diverges. 
Changing variables according to 

— u'^^ 

2/» := 7 > « = 1,...,P, 



(1.8) 



which implies - 

{A)P 

with 




Vi exp 



dy, 



(1.9) 



exp 



1=1 



dy. 



Thus, (|1.9p may be interpreted as the expectation value of a sequence of inde- 
pendent identically distributed (i.i.d.) variables: 



{A)p = E 



where Yi is a sequence of i.i.d. variables distributed according to the law d^{y) = 
exjp{-/3W{y))dy, and Z = J^e^p{—/3W{y))dy. For finite (and large) P, 
the quantity {A)p is not particularly easy to compute. However, its limit when 
P — > oo is easy to compute. Indeed, applying the law of large numbers |35| . 

we know that — Yi converges to the mean E{Yi) of the random variables Yi, 



1=1 



and we hence have (see |6]): 



\im (A) p = A{y*), y*:=E{Y,) = 



yex.p{-l3W{y))dy 
exp{-pW{y))dy 



(1.10) 
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Remark 1 Choosing A{y) — y in the above relation, we observe that y* sat- 
isfies y* = lim (u^)p. In the specific one- dimensional nearest-neighbour case 

discussed here, we actually have y* — {u^)p for any P. Likewise, for any P 

and any 1 < i < P , we have y* = {yi)p = ( )p. 

h 

The property Ijl.lOl) holds under mild assumptions on A and W. The only 
important constraint is that exp{—/3W) should be integrable over the real line. 
Actually, finer results may be easily obtained, as for instance Theorem 1 of [B] . 

In the case of next-to-nearest-neighbour interaction, it is possible to apply 
the same kind of strategy, but it is much more involved. Assume now that, 
instead of (|1.7p . the energy reads 



p-1 



p-1 



(1.11) 



Again, we assume that the observable A depends only on the right-end atom 
position , and that u''^ — 0. Then, using the same change of variables (|1.8|) . 
we have 



1 
1 



A (u^) exp [~[3E{u)] du 



A 



, 1=1 



Vi exp 



X exp 



P-1 

^E^2(2/. + z/,+i; 



i=i 



dy. (1.12) 



Here, instead of i.i.d. variables, we may identify a Markov chain structure. 
Introduce indeed 



k{a, b) — exp 
Then we can recast (|1.12[) as 



■^Wi{a)-^Wiib)-(3W2ia + b) 



:Wi{yi 



/3 



Wiiyf 



p-1 



l[k{y,,y,+,)dy. (1.13) 



i=l 



For the sake of simplicity, assume now that J k{a, b) db = 1. Then we can intro- 
duce a Markov chain {i^i}j>o of kernel fc, namely such that k{a, b) = P(Ki+i = 
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b \ Yi = a), and recast (I1.13P as 



E 




exp (-|vKi(ri)-fiyi(yp)) 


E 


exp 


|w^i(n)-|w^i(rp)) 





Applying next the law of large numbers (this time for Markov chains, rather 
than for i.i.d. variables as previously), we can identify lim {A)p. 

P— s-oo 

In general, the quantity / fc(a, b) db depends on a, and is not equal to 1. So 

the above discussion needs to be amended. It turns out that it is nevertheless 
possible to identify a Markov chain structure in (|1.12p . Hence, applying asymp- 
totic theorems (namely law of large numbers type results) for this setting (see 
e.g. [35]), we have (see [6] for the details): 



hm {A)p = A{y*), 

P— S-CX) 



y 



2/('0i(y)) dy, 



(1.14) 



where -0^ is the invariant measure of the transition kernel of the Markov chain, 
that is, the solution of 



Vipi — Xipi, A — max Spectrum('P), 
where the transition operator V is defined on functions </> : 



(1.15) 



by 



(P^) (x) = / 0(2/) exp 



4wiix)-^Wi{y)~pW2{x + y) 



dy. 



Note that, using standard tools of spectral theory, it is possible to prove that, 
under some reasonable assumptions, problem (|1.15p has a unique solution (we 
refer to [B] for the details). 

Remark 2 As in the nearest-neighbour case, choosing A{y) — y in (|1.14p . we 
have y* — lim (u^)p. Note that we also have 

P-i-oo 



y* = lim 



,1-1 



lim 



for any 1 <^ i <C P. Otherwise stated, y* is the average (reseated) elongation of 
any bond which is in the bulk of the chain. 

Before moving on to the two-dimensional setting, let us briefly describe one 
application of the above strategy. We have considered until now chains of atoms 
at rest, that is, submitted to no external force. The energy of the system 
reads PTT)) (or pTTT]) ). Choosing A{u) ^ u in ([TTTU|) (or pTTil) ). we see that y* is 
the average elongation of the system, at a given finite temperature (macroscopic 
equilibrium length) . It turns out that the above strategy can be easily extended 
to treat the case when some external force is applied to the right-end atom of 
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the chain (at position u^), while we again impose the position of the left-end 
atom: = 0. In this case, the energy, rather than (|1.7I) . reads 

p-i 

i=0 

and the average length of the chain, at the inverse temperature /?, and when a 
force / is appUed at the right-end, reads 

(A^p = ^ [ A (uP) exp [-/3Ef{u)] du, 

for A{u) = u. Following the same arguments as above (see for details), we 
can compute y*c = lim {u'^)p, which is the average elongation of the system 

when submitted to the external force f. We thus have computed the constitutive 
law of the chain (namely the relationship between force and elongation), at finite 
temperature. 

1.2 Setting of the problem 

We now turn to describing the problem we address in this article. The above 
treatment of the one-dimensional setting will be used as a bottom-line to treat 
the case at hand here, namely two-dimensional membrane-like models. In con- 
trast to the one-dimensional case, we are not able to prove the convergences 
we claim. However, we believe that the one-dimensional proofs are sufficient to 
allow for the approximations we are going to make. Extensive numerical results 
reported in Section [3] below show the accuracy of our approximations. 

We consider a two-dimensional atomistic system, where each atom has a 
scalar degree of freedom. The mechanical interpretation is that the system is a 
membrane, and the degree of freedom is the height of the atom (see Figure [1} . 

We hence consider a two-dimensional set of P = (A^ + 1)^ atoms at positions 

u^'^eM, 0<i,j<N. 
We assume a nearest-neighbour interaction, and thus write the energy as 

i<fc<i+l, i<'<j + l 

Note that, as in (ll.7p . we have rescaled the relative position by the factor h = 
1/N. Likewise, we consider in the above sum oriented bonds, from left to right or 
bottom to top (but not the contrary). Finally, as in the one-dimensional case, we 
do not have to assume W to be an even function. The energy (|1.16p is translation 
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Figure 1: Schematic representation of the system: atoms are set on a two- 
dimensional lattice, and have a unique degree of freedom, which is their height. 



invariant, so we fix in the sequel = to remove this invariance. Our aim 
is to compute the canonical average of an observable A that only depends on 
u^'^ , with respect to the Gibbs measure associated to the energy (|1.16p . This 
quantity reads 

{A)m = Z^^I A{u''''')tIn{{u''']) n (1.17) 



where 



JiM{{n'^'])=^M-PE{{u^^^])) (1.18) 
is the non-normalized Gibbs measure, and 



0<i,j<N, (ij)#(0.0) 



is the normalization constant. Our exact aim is to compute {A) ^ in the ther- 
modynamics limit, namely lim {A)^- 

As in [6 , we aim at describing a numerical strategy to compute lim (^4) 

AT— )-oo 

that builds on the specificity that the observable A under study depends only 
on u^'^ , and not on all the atom positions {u*'^ j<N' ^® ^^^^ below, 
we are also able to handle observables that depend only of u'~''^ or u^'", and, 
more generally, of u"'^, u^'" and u^'^ . We build this strategy in Section^ and 



9 



obtain the algorithm summarized in Section [331 Numerical results are reported 
in Section [3l 

Remark 3 In dimension one, in addition to computing lim {A)p (as recalled 

P— i-oo 

in Section ] it is also possible, using a large deviation principle strategy, to 
compute the coarse-grained energy Ecg defined by (jl.6l) . We refer to 'Jy, 2^ for 
the details. In dimension two, such a problem is known to be very difficult, and 
very few theoretical results are known. See for instance [HI [2^ \3^. This 
is why we concentrate here on the simpler question of computing the average of 
an observable. 



2 Numerical strategy 

This section is dedicated to the construction of our numerical strategy. It is 
quite involved from the technical viewpoint. For the sake of clarity, we have 
summarized in Section 12.31 below the resulting algorithm which has eventually 
to be implemented. 

We first rewrite the energy ()1.16|) using the increments between the atoms. 
Introducing 



h 



1 <i <N, 0<j<N, 



= , 0<i<7V, l<j<iV, 

we have 

JV AT 

E ^K.) + w^(?^^j)+E^(^^.o)+E^(yo.j)- (2.1) 

l<i,j<N i=l j=l 

This change of variables is not one-to-one: {N + 1)^ — 1 variables are involved 
in (|1.16p . whereas 2iV^ -I- 2N variables are involved in (|2.1[) . The difference is 
accounted for by imposing the geometric constraints 

for 1 < i,j < N. 

In the variables (xij, yij), the non- normalized Gibbs measure Jl^ reads 

X n n (2.2) 

l<i<Ar i<j<N 

with 

fix) ^eM- mix)), 
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and we recast (|1.17l) as 



V 1=1 / l<i<Af i<i<JV 

X II /(^^j) /(^'j) '^0 + yi-i,j - yi,j - Xij-i) . (2.3) 

l<i,j<N 

Based on our one-dimensional study (recalled in Section [l.ip . our first assump- 
tion is 

[HI] the boundary terms of the Gibbs measure ,„ .s 

(second line of (12. 2p ) are disregarded. 

Thus, we approximate (|2.3p by 

{A)n « (A)N°Bou„dr^ (2.5) 

where 



/ /t\NoBoundr 

with 



{{xij} , {Vi,]}) = W f{xi,j)f{yi,j)So {xij + yi-ij - yij - Xij-i) . 

l<i,j<N 

Again based on the one-dimensional study, we assume that we can apply (|1.14p 
and Remark [21 so that 

[H2] lim {A)fj°^°'^'"^'-' = A{x* +y*), (2.6) 

where 

X* = lim Z]~;^ Xk,i UN {{xij} ,{yi.j}) , (2.7) 

A*— >+oo J 

where k and I are any integers such that 1 <^ fc, Z <C iV. We implicitly assume 
in ([2?7l) - ([2?8| that x* and y* do not depend on fc, provided 1 <C fc, ? < iV. 
Collecting (j2.5p and (j2.6p . we thus write 

lim (^)a, « A(a;* + y*). (2.9) 

All the sequel of this section is dedicated to the computation of (an approxima- 
tion of) X* and y*. 
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Remark 4 Based on the knowledge of x* and y* , we can compute A{x* +y*), 
which is assumed to be an approximation of the average of A{u^'^), according 
to (j2.9p . Using the same assumptions, our strategy leads us to approximate the 
average of A(u^'^) and A{u^'^) hy A{x*) and A{y*) , respectively. We will show 
in Section\3l that A{x* + y*), A{x*) and A[y*) are indeed good approximations 
of the average of A{u^'^), A{u^''^) and A{u'^'^) , respectively. 

As pointed out above, our strategy can straighforwardly be extended to the 
case when the observable A depends on u^'^ , u^'" and u^'^ , although we have 
not numerically tested such cases in the sequel. 



2.1 Numerical computation of x* and y* 

Our strategy consists in setting an appropriate Markov chain structure, as we 
did for the one-dimensional next-to-nearest-neighbour interaction case (see Sec- 
tion [TTT]). The variables Xij in (|2.7p are indexed by {i,j), which belongs to a 
subset of Z^, on which there is no natural ordering. Rather than considering 
the variables Xij, we are going to consider the row of variables {3^i,j}i<i<7V' 
indexed by j. The advantage is that there is a natural ordering between these 
rows of variables. 

We now proceed in details. For any 1 < j < N , define 

We thus have 

TV 

/XAT = n F{x.,^j) F{y,j) 6oix..^j + Ty,j - y,j - a;: j-i), (2.10) 
i=i 

where r is the Bernoulli shift, i.e. (ry)i — j/i-i, and F is the direct product 

JV JV 



i=l 



With this notation, (12.71) reads 



N 

^* " / ^'^^^ TlPi^-:]) Piy-;]) + Ty:,j -y-.j -X..j-l). (2.11) 

We now follow a strategy based on the so-called transfer operators (or transfer 
matrices; see for instance [321 Chapter I, Section 12]). Let us define the transfer 
operator V by 



{'Pip){x,y) = F(x)F{y) ip{t, z) So{t + tz - z - x) dz dt 

= F{x)F{y) / (p{x + z — TZ, z) dz 



12 



for any function 93 defined on x R^. The adjoint operator reads 
{V*'p){x, y) = F{x + Ty~y) / F{z) tp{x + Ty-y, z) dz. 



We can compute the law of {x;^i,y;j) by integrating (|2.10p over aU variables 
except (x-^i^y-^i). The interest of introducing V and V* is that we can now 
recast (|2.1ip in the simple form 

Xk,l Hl{X;^l,y;j) dX;^l dy ; 

m{x:,uy:,i) dx.^i dy.,j 



x*^ lim ^^^^ (2.12) 



with 

MX:,l,y:,l) = (T'^-Vo) iX:,l,y:,l) {{Vn'^'vi) : .U V : .l) , (2.13) 

where ^o{x:^N,y:,N) — F{x;^n) F{y-^pf) and 931 — 1. Assuming N large and 
1 <^ I <^ N, the iterations of V and P* converge, up to renormalization, to the 
spectral projection on the highest eigenvector of V and V* , respectively. Hence, 
the law of {x-j^y-^i) is 

t N ^J'{x■.J,y■.,l) n*{x..,i,y..j) 

MX:.,hy:,l) ^ J 7 ^ \-j (2.14) 

where n and fi* are functions defined on X R^ satisfying 

^fT;^^' A supSpectrum(P) = supSpectrum(7'*). (2.15) 

Remark 5 The operators V and V* depend on N , since they act on functions 
of2N variables. Hence, the highest eigenvalue A as well as the associated eigen- 
vectors (1 and /i* depend on N . To simplify the notation, we have not made this 
dependence explicit in (j2.15p . 

Assuming that the highest eigenvalue of V is simple and isolated, the equa- 
tion ((2?T2]) yields 



lim 



Xk.i n{x.,,i,y.,j) n*{x;,i,y;j) dx.,,i dy-.j 



KX:,hy:,l) fJ-*iX:,hy:,l) dx.,^l dy.,^l 

"xR" 



lim ^-^^ . (2.16) 



Xk t^{x,y) H*{x,y) dx dy 
l^{x,y) y) dx dy 

In the above integrals, x denotes the vector {xi, . . . ,xn) G R^, and 1 ^ k ^ N 
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Remark 6 Note that approximating (|2.12p - (|2.13p hy (j2.16p consists in passing 
to the limit N oo faster in the vertical direction than in the horizontal direc- 
tion. Indeed, the operator V acts on functions of 2N variables, which correspond 
to increments {xi.j,yi.j}^^J^ ^, for some j. These are all the increments along 
a row of atoms in the horizontal direction. We hence freeze the number of par- 
ticles in the horizontal direction and pass to the limit in the vertical direction, 
since we consider the eigenmode of highest eigenvalue of V and V* . 

In order to take benefit from (and proceed furtlier tlian) (|2.16p . we need to 
solve the eigenvalue probiem (|2.15p . This problem is posed on functions on 2N 
variables. Solving this problem directly is hence out of reach numerically. The 
approach we suggest consists in further approximating the problem, arguing 
somewhat similarly to a mean-field theory. Using a sequence of simplifying 
assumptions (see in particular (j2.19l) below) , we are going to derive an integrated 
form of (I2T5)) (see (|Z!20)) and ((2?2T]) below). Eventually, this will allow us to 
recover a Markov structure on low dimensional objects. The invariant measure of 
such Markov chain will hence be the solution of an eigenvalue problem, similarly 
to (I2.15p . but, in contrast to (I2.15p . in a low-dimensional setting, and thus 
amenable to numerical computations. 

Before proceeding along these lines, we make the following observation. We 
note that any /i satisfying (|2.15p is of the form 

V(a;, y) e X R^, ^{x, y) = F{x)F{y)v{x), with 

v{x) = / F{z)F{x + z — Tz)i'{x + z ~ Tz)dz. (2.17) 

Similarly, any /i* solution to (I2.15P is of the form 

V(a;,2/) G R^ X R^, n* {x,y) ^ F{x + ry - y) v* [x + ry - y), with 

VxGR^, \v*{x)^ F{z) F{x + TZ - z) v*{x + TZ - z) dz. (2.18) 

In the sequel, we describe a strategy to approximate v and v* , which is based 
on a direct product ansatz. The functions v and v* depend on fewer variables 
than and /^*, and are hence easier to handle. 

2.2 A direct product ansatz 

Considering that the problems (I2.17P and (I2.18P are (at least asymptotically) 
symmetric in the variables, we introduce the ansatz 

N N 

[H3] v{x) ^\{g{x,), v''{x)^\{g*{x,), (2.19) 

i=l i=l 
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insert it in (|2.17p and (|2.18p . and integrate with respect to all Xi except one. 
This gives the following eigenproblems on g and g*: 



^g{xi) = / f{.Xi + Zi- Zi^i)f{zi)f{zi^i)g{xi + Zi- Zi^i)dzi^idzi, (2.20) 
Xg*{xi) = / f{xi + - Zi)f{zi)f{zi^i)g*{xi + - Zi)(izjdzi_i (2.21) 



where A is the largest possible eigenvalue. Note that (I2.20p and (|2.2ip are the 
same equation. Using the Krein-Rutman theorem |32| . it is possible to prove 
that for many types of interactions, the corresponding eigenvalue is simple. 
Hence, g — g* ■ In addition, observe that (I2.20p is an eigenvalue problem on 
functions of one variable: it can be solved numerically. 

Using (|2.17p . (|2.18p and (|2.19p . we thus have, up to a normalization factor, 
that, for any {x,y) S M^^, 

N N 

l^(.x,y) = Ylfixj)fiyj)g{xj), n*{x,y) = J| /(a;^ +yj+i -yj)g(a;j +yj+i -%), 

(2.22) 

where g is the solution of (|2.20p . 

Remark 7 Note that the direct product ansatz (j2.19l) does not imply that, in the 
law (j2.14p of {x;^i,y;^i), the increments {xj_i,yj^i), j ,N, are independent. 

Indeed, in view of (j2.22p . the numerator of (j2.14l) reads 

l^{x:,uy:,i) fi*{x..j,y.,j) = 

N 

which is not a direct product of functions depending on {xjj,yjj)jLi. 
We now insert (|2.22p in (|2.16p and obtain 



^* = ^^^^^ ^Jv / T\f{xj)f{yj)g{xj) 

N-^+QC JrNxE« J-Jl 

X fixj + Vj+i - yj)g{xj + yj+i - dxdy, (2.23) 

N 

where Zn = TT f{xj)f{y.j)g{xj)f{xj +yj+i-y.j)g{xj -y^) dxdy. 

The Markov structure present in the right hand side of (I2.23P now allows us 
to proceed. We introduce the transfer operator Q defined by 



{Q^){xj+i,yj+i) = fixj+i) f{yj+i) / g{xj) 

JR2 

X f{xj + yj+1 ~ yj) g{xj + y^+i - y^) ip{xj,yj) dxj dy^, (2.24) 
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its adjoint Q* , and the invariant measure equations 

{ st/^^ityV, ^ = supSpectrum(Q) = supSpectrum(Q'^). (2.25) 
Using Q and Q*, we recast (|2.23p as 



Xk (Q Vo)(a;fc,yfc) ((Q*) ^i^i) {xk.Vk) dxkdyk 
X* = lim ^ii^^^ (2.26) 



N 



(Q'^Vo) (a;fe,yfe) ((Q*)'^-Vi) {xk,yk) dxkdyk 



for some ^/iq and ipi . Using again the Krein-Rutman theorem, we obtain that the 
highest eigenvalue of Q is simple and isolated, hence the iterations of Q converge, 
up to renormalization, to the projection on ijj, the eigenvector associated to the 
highest eigenvalue. We hence infer from (|2.26p that 



Xk'ipixk,yk)il^*{xk,yk)dxk dyk 
^{xk,Vk) ^*{xk, Vk) dxk dyk 

2 

yk ip{xk, yk) 'ip*{xk,yk) dxk dyk 
^{xk.yk) V{xk,yk) dxk dyk 



(2.27) 



Note that, owing to the specific structure of Q, its eigenvector i(j{x,y) reads 
y{x,y)eRxR, ^{x,y)^ J{x)f{y)^{y), (2.28) 
where, for any y' e M, 

mPiy')^ [ 9ix)f{x + y'-y)g{x + y'-y)fix)f{y)^{y)dxdy. (2.29) 

Hence, computing the function ip, which depends on two scalar variables, amounts 
to determining tp, namely solving an eigenvalue problem on functions of one 
variable. In a similar way, the eigenvector ip*{x, y) of Q* reads 

W{x,y)eRxR, ij*ix,y)=g{x)i^ix-y), (2.30) 
with, for any a; £ R, 

r]ip*{x)= / f{x' + y') f{y') f{x + y')g{x + y')g{x' + y')ip*{x') dx' dy'. (2.31) 



Again, we are left with solving an eigenvalue problem on a function of one 
variable. 
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2.3 Summary 

The above approach results in an algorithm for computing an approximation of 
(A)pf. This algorithm may be summarized as follows: 

1. compute the function g : M h- > R satisfying the eigenvalue equation (j2.20p . 
Note that A should be the largest eigenvalue of the corresponding operator; 

2. with g, compute the operator Q given by (|2.24p . and compute the func- 
tions tp and solution to (|2.25l) ; this amounts in practice to solving the 
eigenvalue problems (P?^ and ((OTt for : IR i-^ K and : R R, 
where again ry is the largest eigenvalue; the functions ip and ip* are next 
computed using ([2?28| and ([230| : 

3. with tp and tp*, compute x* and y* given by ()2.27p : 

4. approximate the limit of {A)^ by A{x* + y*), following (|2.9I) . 

Remark 8 yls pointed out in Remark[^ our strategy can also handle observ- 
ables that depend on u^'^ , u^'^ and m*^' . We then take the approximation 

hm {A(u^^^,u^''',u'''''))m - A{x* + y\x\y*). 

2.4 The Gaussian case 

In this section, we consider the special case when, in (jl.l6p . we choose 

W{x) = x^, fix) = exp(-;3M^(a;)) = exp(-;3a;2). 

Let us show that, in this case, our strategy formally yields the correct value 
for X* and y* . To simplify the notation, we assume henceforth that the inverse 
temperature is pi — 1. 

First, since hn is an even function, we infer from (|2.7p that x* ~ 0. Let 
us now follow our strategy, as summarized in Section 12.31 above. Taking the 
Fourier transform of the equation ()2.20|) . we have 

Xgico) = f{uj)J{-u:)T9{^). 

This is an eigenvalue problem, and we are interested in the case when A is 
maximal. In such a case, and under mild assumptions on / (see 6J for the 
details) , one easily proves that the corresponding solution g is unique and never 
vanishes. In addition, any other eigenvector must vanish at some point of R 
(see O Chap. 1, Theorem 1.2]). Hence, if we find a non- vanishing solution 
of the above problem, it is the unique solution we are looking for. Postulating 
that g{x) = exp(— ax^) for some a, we deduce from the above relation that a 
should satisfy 1 = 2a + 2o? . This equation has two real solutions, and only one 
of them is positive. It thus determines a > 0. This gives us a positive solution, 
which, according to the above argument, corresponds to the largest possible A. 
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We now consider the eigenvalue problem (I2.29p . Taking the Fourier trans- 
form of this equation, and postulating that ipi^) — exp(— tx^) for some t, we 
find that r should satisfy 1 + a = 2t + 2t^. Again, this equation has two real 
solutions, and only one of them is positive. It thus determines r > 0. 

We finally turn to the eigenvalue problem (|2.3ip , and assume that its solution 
reads V (a^) = exp(— 70;^) for some 7. A tedious computation then shows that 
7 should satisfy 

1 1 1 

- = + > 0. 

7 1 + a 1 + T 

Having identified the functions g, "0 and "0 , we can now evaluate the right-hand 
side of (|2.27p . Since /, and are even functions, wc note that 

vWix, y) = f{x)f{y)tp{y) g(x)0*(a; - y) = ^{-x, -y). 

Using the change of variable {x, y) n> (— x, —y) in the right-hand side of (|2.27p . 
we deduce that 

xip{x,y)ip*{x,y) dx dy 

^ = 0. 

-0(a;, y) ip* (x , y) dx dy 

2 

Consequently, our strategy yields the exact result, namely a;"^ = 0. In the 
following section, we numerically consider a more generic case, for which no 
analytical solution is available. 



3 Numerical results 

For our numerical simulations, we choose the interaction potential 

W{x) = ^ix~iy + ^x\ 

Note that W is not an even function, neither does it satisfy / x exp{— (3W (x)) dx — 

0. Hence the difference in height over a given edge of the membrane lattice does 
not average to 0. For a given choice of the inverse temperature /3, we have 
followed the strategy summarized in Section [^751 that yields x* and y* defined 
by (|2.27p . This allows to compute the reduced model results, that is, the quan- 
tity A{x* + y*), for any given observable A. 

The reference model result is {A)n defined by (|1.17p . To compute this 
quantity, we resort to the stochastic differential equation 

dut = -VuE{ut) dt + dBt in m(^+i)'~\ (3.1) 

where Bt is a standard Brownian motion in dimension (N A- 1)^ — 1. Of course, 
this is a very expensive way of computing {A) pf , and we use it only to validate 
our approach, which is far less computationally demanding (see Table [1] below. 
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where we compare computational costs). As recalled in the Introduction, un- 
der mild assumptions on the potential W and the observable A, that are here 
satisfied, an ergodic theorem holds, and the canonical average (|1.17p is given by 

iAU=lunJl"A{ur)dt, 

for (almost) all initial conditions ut=o- In practice, the equation p.ip is numeri- 
cally integrated with the forward Euler scheme (also called the Euler-Maruyama 
scheme) 

um+i = Urn -At V„£;(u„) + V2At^-i G„ (3.2) 

where Gm is a vector of Gaussian random variables in dimension {N -1-1)^ — 1, 
and At is a small time step. In turn, the canonical average is approximated by 

1 

(^)--Ji5LME^m- (3-3) 

m— 1 

In practice, we have integrated p.ip (using the scheme (|3.2p ') for a large but 
finite number M of time steps. We thus take the approximation 

M 

(^)a^«^E^("™'') (3-4) 

rn — l 

for large M and small At. The deterministic quantity (A) at is thus approxi- 
mated by a random number. To compute error bars for {A)[^ (that is, confi- 
dence intervals), we have simulated many independent realizations of the dy- 
namics p.ip . 

The reference model results reported here have been obtained using (|3.2p - 
(|3.3p with At = 10^'^ and a number of time steps M large enough such that 
convergence is reached in p.3p . up to statistical noise. 

3.1 Average of observables: convergence with 

In Figures [2] and m we compare the reference observable average {A)pf with its 
approximation A(x* + y*)., for several observables A, and for increasing values of 
N (the temperature is fixed at /3~^ = 1). We observe a good agreement between 
A{x* + y*) and {A)n, the latter seeming indeed to converge to the former, when 
N ^ (x: for iV = 100, the relative difference is of about 1%. 

In Figure m we compare {u^'^)n with x* , in line with Remark [H We 
again observe a good agreement between the reduced model and the reference 
model results. On the same figure, we compare {u^'^)n with y*, with similar 
conclusions. 

For the sake of completeness, we compare in Table[l]the costs for computing 
the reference value {A) with those associated to the reduced model approach. 
We clearly see that the latter are much smaller than the former. 
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Figure 2: Comparison of (A) n (reference result, defined by ()1.17|) ). with A{x* + 
y*) (reduced model result), for several values oi N [jS = 1; left: A{u) = v?\ right: 
A{u) = u). 
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Figure 3: Comparison of {A)n (reference result defined by (|1.17|) ). with A{x* 
y*) (reduced model result), for several values of N {A{u) ~ exp(u), f3 ~ 1). 





iV 25 


TV = 50 


N = 100 


Reduced model 


CPU time (days) 


0.2 


0.7 


4.6 


0.025 



Table 1: First 3 columns: CPU time for computing the reference value {A)n 
(for one choice of observable A, at one given temperature) using the approxi- 
mation (|3.4p along one single realization of p.2p . with At = 10~^ and M = 10^ 
time steps. Last column: CPU time needed by the algorithm of Section [2T3l (in 
practice, eigenproblems have been solved on the interval / = [—4; 5], using a 
spectral basis of 20 functions; integrals on the interval / have been computed 
using the 6 Gauss points quadrature rule, discretizing / into 20 subintervals). 
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Figure 4: Left: Comparison of {u'^'°)n (reference result) with x* (reduced 
model result), for several values oi N {f3 = 1). Right: Comparison of (w°'^)Ar 
(reference result) with y* (reduced model result). 



3.2 Average of observables: temperature dependence 

In Figures [5] and ini we compare the reference observable average {A)^ with its 
approximation A{x* +y*), for several observables A, and for different tempera- 
tures (the reference system is run with N — 100 particles in each direction). We 
observe a good agreement between A{x*+y*) and {A)i\f, for all the temperatures 
considered. 




1/13 IIP 

Figure 5: Comparison of (A)jv (reference result defined by ()1.17p ). with A{x* + 
y*) (reduced model result), as a function of the inverse temperature /3 [N = 100; 
left: A{u) = u^; right: A{u) = u). 

On Figure [71 we compare {u^'°)n (respectively {u°''^)n) with x* (respec- 
tively y*). We again observe a good agreement between the reduced model and 
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0.25 0.5 0.75 1 



V/3 

Figure 6: Comparison of {A)n (reference result defined by (I1.17p ). with A{x* + 
y*) (reduced model result), as a function of the inverse temperature /3 {A{u) = 
exp(M), N = 100). 

the reference model results, for all the temperatures we have tested. 




1//3 

Figure 7: Left: Comparison of {u^''^)n (reference result) with x* (reduced 
model result), as a function of the inverse temperature /3 {N = 100). Right: 
Comparison of (u°'^)jv (reference result) with y* (reduced model result). 



3.3 Correlations 

In this section, we compare the correlations of successive increments in the 
membrane. More precisely, we compute, according to the reference model and 
the reduced model, the average of products of the type Xk,iXk+p,i for some p > 0. 
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Recall that 

Xk.l = 



h 

is the (rescaled) difference in height between atoms at sites {k — 1,1) and (fc, I). 
This kind of tests are more demanding than the ones considered above. 

The reference model computations have been performed on a system of {N + 
1)^ atoms. The correlation can be defined from all the increments except those 
at the boundaries, or just the increments at the center of the system. Both 
yield very close results (the difference is smaller than 0.2 %). In addition, those 
reference computations have been computed by averaging several independent 
realizations. The number of realizations we have considered is large enough so 
that the statistical noise is smaller (in relative error) than 0.2 %. Results given 
by the reference computations are gathered in the second column of Table [21 

The reduced model values have been computed by starting from the exact 
values, e.g. 

(4,/> = ^7V^ J , {Vij}) , 

and following the steps described in Section This yields the quantity 



X,, ii(xk,yk) ip*{xk,yk) dxk dyk 

(^M)rcd = '-^r (3-5) 

ipix k,yk) ip*ixk,yk) dxk dyk 



as an approximation of {x^. In turn, the quantity {xk i Xk+i i) is approximated 

by 

{Xk,l Xk+iu)rcd — 

XkXk+1 Tpixk^Hk) c{xk,yk,xk+i,yk+i) 4>* {xk+i,yk+i) dxk dyk dxk+i dyk+i 
/ il^{xk,yk) c(xfe , yk , Xfc+i , yk+i) V'* (xk+i, 2/fc+i ) dxk dyk dxk+i dyk+i 

(3.6) 

with 

c{xk,yk,Xk+i,yk+i) 

g(xk) .f(xk+i) fiyk+i)g{xk + yk+i - yk) f{xk + Vk+i - yk)- 

Note that formulas (|3.5p and p.6p are consistent with each other. Indeed, if 
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Xk+i,i is replaced by Xk,i in (13.61) . we have 

{Xk,l 2;fc,;)rcd 

xl ip{xk,yk) c{xk,yk,Xk+i,yk+i)iJ*{xk+i,yk+i) dx^ dyu dxk+i dyu+i 

ip{xk, yk) c{xk,yk, Xk+i,yk+i) ■ip*ixk+i,yk+i)dxk dyt dxk+i dyt+i 
xli>{xk,yk) {Q*V){xk,yk)dxk dyk 

ip{xk,yk) {Q*ip*){xk,yk)dxkdyk 
xl Tpixk^yk) ip*{xk,yk) dxk dyk 
i^{xk, yk) i^*{xk,yk) dxk dyk 
and we recover (13.51). We have used in the above derivation the fact that 



c{xk,yk,xk+i,yk+i)i^*{.xk+i,yk+i)dxk+i dyk+i = 

Q*ip*{xk,yk) = ?7i/'*(a;fe, j/fc). 



Integrals in dimension 2 and 4 need to be evaluated to compute p. 61) (and 
similarly, integrals in dimension 6 need to be evaluated to compute {xk,i Xk+2,i)i-cd)- 
These integrals have been computed almost exactly using numerical quadra- 
tures. We gather the obtained results in the third column of Table [2j 





Reference model 


Reduced model 


Relative difference 




0.50668 


0.52808 


4.05 X IQ-^ 


ivii) 


0.50684 


0.51217 


1.04 X 10"^ 


{Xk,l Xk+ld) 


0.28097 


0.31697 


0.114 


{Xk,l Xk+2j) 


0.31215 


0.33089 


5.66 X 10-^ 


{Xk,l) 


0.57310 


0.581 


1.36 X 10"^ 



Table 2: Comparison of the correlations of successive increments according to 
the reduced model with the reference model values (inverse temperature /? = 1, 
N = 100, 1 <^ k,l N). The reference values are expectations (e.g. of x^, ;) 
with respect to the Gibbs measure Z^^Jl^ defined by (|1.18p . 

We observe that the difference between the reference result and the reduced 
model results is of the order of 5% for all values except one. We consider 
that such an accuracy is good, given the number of assumptions on which the 
numerical strategy leading to the reduced model values is based. 
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